rm(list = ls())
require(zTree)
require(stargazer)
require(mfx)
require(stats4)
require(nlWaldTest)
require("multiwayvcov")
require(msm)
require(miceadds)
require(restriktor)


zTT1= zTreeTables("C:/Users/naten/OneDrive - University of Tennessee/Dean Stuff/Ambueletalexperiments/Data/August2Expansionpath2021.xls")
subjects1=zTT1$subjects

zTT2= zTreeTables("C:/Users/naten/OneDrive - University of Tennessee/Dean Stuff/Ambueletalexperiments/Data/September21ExpansionPath2021.xls")
subjects2=zTT2$subjects
subjects2$Subject=subjects2$Subject+100

zTT3= zTreeTables("C:/Users/naten/OneDrive - University of Tennessee/Dean Stuff/Ambueletalexperiments/Data/September 29ExpansionPath2021.xls")
subjects3=zTT3$subjects
subjects3$Subject=subjects3$Subject+200

zTT4= zTreeTables("C:/Users/naten/OneDrive - University of Tennessee/Dean Stuff/Ambueletalexperiments/Data/October13expansionpath2021.xls")
subjects4=zTT4$subjects
subjects4$Subject=subjects4$Subject+300


subjects=rbind(subjects1,subjects2,subjects3,subjects4)

#Clean data
subjects$responder=(subjects$"response1[50]">0)
subjects=subset(subjects,responder>0)

uidlist=unique(subjects$Subject)
incentivec=c(5,95)

#create a frame with observations
obsframe=data.frame(uid=vector('numeric'),block=vector('numeric'), incentive=vector('numeric'), number=vector('numeric'), state=vector('numeric'), response=vector('numeric'),react=vector('numeric'))
for(i in 1:length(uidlist)){
playframe=subjects[i,]
uidplace=uidlist[i]

#block 1
for (j in 1:50){
blockplace=1
incentiveplace=incentivec[playframe[1,683]]
numberplace=j
stateplace=playframe[1,(736+j)]
responseplace=playframe[1,(72+j)]
reactplace=playframe[1,(122+j)]
addframe=data.frame(uid=uidplace,block=blockplace, incentive=incentiveplace, number=numberplace, state=stateplace, response=responseplace,react=reactplace)
obsframe=rbind(obsframe,addframe)
}

#block 2
for (j in 1:50){i
blockplace=2
incentiveplace=incentivec[playframe[1,684]]
numberplace=j
stateplace=playframe[1,(786+j)]
responseplace=playframe[1,(222+j)]
reactplace=playframe[1,(272+j)]
addframe=data.frame(uid=uidplace,block=blockplace, incentive=incentiveplace, number=numberplace, state=stateplace, response=responseplace,react=reactplace)
obsframe=rbind(obsframe,addframe)

}

}



#define correctness
obsframe$correct=(obsframe$state==obsframe$response)

#mean correctness
correct5=sum(obsframe$correct*(obsframe$incentive==5))/sum(obsframe$incentive==5)
correct95=sum(obsframe$correct*(obsframe$incentive==95))/sum(obsframe$incentive==95)
correctness=c(correct5,correct95)

obsframe5=subset(obsframe,incentive==5)
obsframe95=subset(obsframe,incentive==95)

#change terminology to work with older scripts
expframe=obsframe

expframe$incentive5=(expframe$incentive==5)
expframe$incentive95=(expframe$incentive==95)

expframe$Bresponse <-expframe$response==2
expframe$Bstate <- expframe$state==2
expframe$Aresponse <-expframe$response==1
expframe$Astate <- expframe$state==1

expframe5=subset(expframe,incentive==5)
expframe95=subset(expframe,incentive==95)


#write.table(expframe, "C:/Users/naten/OneDrive - University of Tennessee/Dean Stuff/Revision/Reports/AmbhuelExpansionPath.csv", sep = ",", col.names = T, append = F)

#NIAS in aggregate

#NIAS in aggregate
#Base Regression
expframe$incentivedec=expframe$incentive/100
modelNIAS=lm(Aresponse~factor(Astate+incentivedec)+0,data=expframe)

#Adjust COV
vcovmodNIAScluster=cluster.vcov(modelNIAS, expframe$uid)

modelniassterr=vector("numeric")
#Model Nias sterrors
for(i in 1:4){
modelniassterr[i]=sqrt(vcovmodNIAScluster[i,i])
}

logitNIAS=glm.cluster(data=expframe,Aresponse~factor(Astate+incentivedec)+0, cluster="uid", weights=NULL, subset=NULL,  family="binomial" )

#Check NIAS and NIAC on the individual level
#NIAS


NIASAframe=data.frame(uid=vector(),A_state_avg_effect=vector(),Pvalue=vector(),incentive=vector())
for (i in 1:length(uidlist)){
for(j in 1:2){
addframe=0
playframe=subset(expframe,uid==uidlist[i] & incentive==incentivec[j])
if(sum(playframe$Aresponse)==0|sum(playframe$Bresponse)==0){
addframe=data.frame(uid=uidlist[i],A_state_coefficient=0,Pvalue=1,incentive=incentivec[j])
NIASAframe=rbind(NIASAframe,addframe)
}else{
model=logitmfx(Aresponse~Astate,data=playframe,robust=TRUE) 
addframe=data.frame(uid=uidlist[i],A_state_coefficient=model[[1]][1],Pvalue=model[[1]][4],incentive=incentivec[j])
NIASAframe=rbind(NIASAframe,addframe)
}
}
}

#Check NIAC 

restrictmat=matrix(c(0,1),nrow=1,ncol=2, byrow=TRUE)

NIACindiframe=data.frame(uid=vector('numeric'),fstat=vector('numeric'),pvalue=vector('numeric'))
for(k in 1:length(uidlist)){
playframe=subset(expframe,uid==uidlist[k])

playmodel=glm(correct~incentive95,data=playframe,family=binomial())
playtest=conTest(playmodel,constraints=restrictmat,type="B",rhs=c(0))
addframe=data.frame(uid=uidlist[k],fstat=playtest$Ts,pvalue=playtest$pvalue)

NIACindiframe=rbind(NIACindiframe,addframe)
}
sum(NIACindiframe$pvalue<0.05)/length(NIACindiframe$pvalue)

#NIAS and NIAC together
comboframe=data.frame(uid=vector('numeric'),NIASviol=vector('numeric'),NIACviol=vector('numeric'))
for(i in 1:length(uidlist)){
NIACplayframe=subset(NIACindiframe,uid==uidlist[i])
NIASplayframe=subset( NIASAframe,uid==uidlist[i])
addframe=data.frame(uid=uidlist[i],NIASviol=(sum((NIASplayframe$A_state_coefficient<0)*(NIASplayframe$Pvalue<0.05))>=1),NIACviol=(NIACplayframe$pvalue<0.05))
comboframe=rbind(comboframe,addframe)
}

comboframeclean=comboframe
NIASonly=100*sum((comboframeclean$NIASviol==1)*(comboframeclean$NIACviol==FALSE))/length(uidlist)
NIAConly=100*sum((comboframeclean$NIASviol==0)*(comboframeclean$NIACviol==TRUE))/length(uidlist)
both=100*sum((comboframeclean$NIASviol==1)*(comboframeclean$NIACviol==TRUE))/length(uidlist)
neither=100*sum((comboframeclean$NIASviol==0)*(comboframeclean$NIACviol==FALSE))/length(uidlist)
violatereport=data.frame(violate=c("NIAS only","NIAC only","Both","Neither"),percent=c(NIASonly,NIAConly,both,neither))

#************************************
stargazer(violatereport,summary=FALSE,rownames=FALSE)
#************************************

goodlist=comboframe$uid*(1-comboframe$NIASviol)*(1-comboframe$NIACviol)

#data cleaning
chooseB=vector('numeric')
for(i in 1:length(uidlist)){
playframe=subset(expframe,uid==uidlist[i])
chooseB[i]=(sum(playframe$Bresponse)>5)
}

zeroedlist=chooseB*uidlist
cleanlist=zeroedlist[zeroedlist>0]






#data cleaning
chooseB=vector('numeric')
for(i in 1:length(uidlist)){
playframe=subset(expframe,uid==uidlist[i])
chooseB[i]=(sum(playframe$Bresponse)>5)
}

zeroedlist=chooseB*uidlist
cleanlist=zeroedlist[zeroedlist>0]



#try with bootstrapping
N=1000
n=length(cleanlist)
randomdraw=sample(1:n,n*N,replace=TRUE)
consframestateA=data.frame(acc5=vector('numeric'),acc40=vector('numeric'),acc70=vector('numeric'),acc95=vector('numeric'))
consframestateB=data.frame(acc5=vector('numeric'),acc40=vector('numeric'),acc70=vector('numeric'),acc95=vector('numeric'))

for(k in 1:N){
bootframe=data.frame()
for(l in 1:n){
addframe=subset(expframe,uid==cleanlist[randomdraw[(k-1)*n+l]])
bootframe=rbind(bootframe,addframe)
}



#for state A
frameA=subset(bootframe,Astate==TRUE)

A1=sum(frameA$correct*(frameA$incentive==5))/sum(frameA$incentive==5)
A2=sum(frameA$correct*(frameA$incentive==95))/sum(frameA$incentive==95)

#for state B
frameB=subset(bootframe,Bstate==TRUE)
B1=sum(frameB$correct*(frameB$incentive==5))/sum(frameB$incentive==5)
B2=sum(frameB$correct*(frameB$incentive==95))/sum(frameB$incentive==95)

addframeA=data.frame(acc5=A1,acc95=A2)
addframeB=data.frame(acc5=B1,acc95=B2)

consframestateA=rbind(consframestateA,addframeA)
consframestateB=rbind(consframestateB,addframeB)
}

#for the lambda estimates

accavgA=c(mean(consframestateA$acc5),mean(consframestateA$acc95))
acchighA=c(sort(consframestateA$acc5)[950],sort(consframestateA$acc95)[950])
acclowA=c(sort(consframestateA$acc5)[50],sort(consframestateA$acc95)[50])

accavgB=c(mean(consframestateB$acc5),mean(consframestateB$acc95))
acchighB=c(sort(consframestateB$acc5)[950],sort(consframestateB$acc95)[950])
acclowB=c(sort(consframestateB$acc5)[50],sort(consframestateB$acc95)[50])

lambdaA=vector('numeric')
lambdasterrA=vector('numeric')
lambdahighA=vector('numeric')
lambdalowA=vector('numeric')

lambdaA[1]=min(1000,0.5*5/(log(accavgA[1]/(1-accavgA[1]))))
lambdaA[2]=min(1000,0.5*95/(log((accavgA[2])/(1-(accavgA[2])))))

lambdahighA[1]=min(1000,0.5*5/(log(acchighA[1]/(1-acchighA[1]))))
lambdahighA[2]=min(1000,0.5*95/(log((acchighA[2])/(1-(acchighA[2])))))

lambdalowA[1]=min(1000,0.5*5/(log(acclowA[1]/(1-acclowA[1]))))
lambdalowA[2]=min(1000,0.5*95/(log((acclowA[2])/(1-(acclowA[2])))))

#for the lambda estimates
lambdaB=vector('numeric')
lambdahighB=vector('numeric')
lambdalowB=vector('numeric')

lambdaB[1]=min(1000,0.5*5/(log(accavgB[1]/(1-accavgB[1]))))
lambdaB[2]=min(1000,0.5*95/(log((accavgB[2])/(1-(accavgB[2])))))

lambdahighB[1]=min(1000,0.5*5/(log(acchighB[1]/(1-acchighB[1]))))
lambdahighB[2]=min(1000,0.5*95/(log((acchighB[2])/(1-(acchighB[2])))))

lambdalowB[1]=min(1000,0.5*5/(log(acclowB[1]/(1-acclowB[1]))))
lambdalowB[2]=min(1000,0.5*95/(log((acclowB[2])/(1-(acclowB[2])))))


lambdas <- rbind(lambdaA,lambdaB)
lambdahigh <- rbind(lambdahighA,lambdahighB)
lambdalow <- rbind(lambdalowA,lambdalowB)


#***********************

barCenters <- barplot(lambdas,
                 names.arg = incentivec,
                  beside = TRUE, las = 2,
                  ylim = c(0, 50),
                  cex.names = 0.75, xaxt = "n",
                  main = "",
			col=c("blue","orange"),
			ylab="Cost",
			xlab="Incentive",
                  border = "black", axes = TRUE,xpd = FALSE, legend=c("state 1","state 2"))
axis(side=1,at=(barCenters[1,]+barCenters[2,])/2,labels=incentivec,  xlab = "Incentive")
box(bty="l")

segments(barCenters, lambdalow, barCenters,
       lambdahigh, lwd = 1.5)

arrows(barCenters, lambdalow, barCenters,
      lambdahigh, lwd = 1.5, angle = 90,
       code = 3, length = 0.05)
#**********************************